In this example, we track the numbers of vehicles on different links in a transportation network from the observations about probe vehicles. We have synthesized the following dataset from MATSIM equil example:
The transportation network looks like the following. Individuals in this scenario leave their homes on link 1 to their workspace on link 20 around 6:30am, and leave their workspace to their home around 7pm. Our goal is to track the numbers of vehicles at different locations in a transportation network from the numbers of probe vehicles at link 1 and link 20 using a particle filter.
synthesized data from MATSIM equil
The particle filter algorithm takes 3 pieces of information as input: - h(x,t) the event rates of the uniformized stochastic kinetic process, - S.uniformized the stoichiometric matrix, and - obs.prob(x,t) the likelihood of latent state with respect to observations.
It generates as output: - Xt_est the estimated numbers of vehicles at different locations.
while(t < nrow(Xt_est)) {
gamma = max(apply(Xt, 1, function(x) {sum(h(x,t))*gamma.inc} ))
event.time.sched <- t+rexp(1,rate = gamma)
# resampling/updating
if ((trunc(event.time.sched)-trunc(t)) > 0) {
for (i in c(1:(trunc(event.time.sched)-trunc(t)))) {
W = obs.prob(Xt,t,i)
Xt = Xt[sample(length(W),particles,replace=TRUE,prob=W),]
Xt_est[trunc(t)+i,] = colMeans(Xt)
}
}
# sampling/prediction
Xt = t(apply(Xt, 1, function(x) {
event.rate = c(h(x,t)/gamma, 1-sum(h(x,t))/gamma)
vt <- sample(c(1:length(event.rate)),1,prob = event.rate)
x <- x + S.uniformized[,vt]
}))
# updating time
t = event.time.sched
}
# plot result
layout(matrix(1:ncol(Xt_real),ncol=1), heights=apply(Xt_real,2,max))
par(mar=c(0,0,0,0),oma=c(5,4,4,2)+.1)
for(i in 1:ncol(Xt_real)){
plot(Xt_real[,i],type='l',col=1,xaxt='n',xlab='',ylab=colnames(Xt_real)[i],lty=2)
lines(Xt_est[,i], col=2,lty=1)
if(i==1) legend('bottomright',col=1:2,lty=2:1,legend=c('groundtruth','particle filter') )
}
axis(side=1)